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Chapter  1 

Introduction 


Under  this  MRI  funded  by  HEL  JTO  and  AFOSR,  UCLA,  in  collaboration  with  Michigan  Tech, 
Georgia  Tech,  MZA  Associates  Corporation,  Tempest  Technologies,  Trex  Enterprises  Corporation, 
ATK  Mission  Research  and  the  Air  Force  Institute  of  Technology  (AFIT),  has  established  a  com¬ 
prehensive  research  program  in  high-performance  control  of  high  energy  lasers  (HEL),  modeling 
and  simulation  of  HEL,  wavefront  sensing,  and  target  tracking.  Under  this  program,  researchers 
have  developed  adaptive  filtering  and  control  methods  for  wavefront  prediction  and  correction  and 
precise  pointing  of  laser  beams  to  compensate  for  the  effects  of  atmospheric  turbulence,  platform 
vibration,  target  motion  and  sensor  noise,  all  of  which  degrade  the  performance  of  laser  weapons 
and  communication  systems. 

Recent  improvements  in  laser  power  and  wave  front  control  technology  for  space  surveillance 
and  laser  anti-satellite  and  anti-ballistic  missile  weapons  has  motivated  interest  in  extensions  and 
alternative  uses  of  this  technology.  Of  particular  interest  are  directed  energy  weapons,  such  as  lasers. 
The  agility  and  speed  with  which  laser  weapons  can  operate,  combined  with  potential  pinpoint  accu¬ 
racy  and  low  collateral  damage  associated  with  these  weapons  make  laser  weapons  highly  desirable 
for  a  variety  of  applications,  including  high  altitude  Airborne  Laser  (ABL),  low  altitude  tactical 
battlefield  scenarios,  and  marine  scenarios.  However,  considerable  fundamental  scientific  work  must 
be  conducted  to  bring  these  weapons  to  the  battlefield  with  the  capability  to  deliver  energy  to  the 
target  through  the  atmosphere  in  each  scenario  of  interest.  Overcoming  the  technological  barriers 
requires  a  multidisciplinary  approach  that  combines  research  on  the  physics  of  beam  propagation 
and  control  as  well  as  on  the  mathematics  an  practical  implementation  of  advance  methods  for  active 
beam  control  and  target  tracking. 

The  research  team  for  this  multidisciplinary  research  initiative  is  making  a  comprehensive,  in¬ 
tegrated  attack  on  the  broad  range  of  modeling  and  simulation,  beam  control,  and  target  tracking 
problems  that  must  be  solved  to  achieve  the  potential  of  high  energy  laser  systems.  This  report 
discusses  the  main  results  of  the  research. 

Members  of  our  team  have  worked  in  close  collaboration  so  that  each  group  can  take  advantage 
of  the  expertise  in  the  other  groups.  The  modeling  and  simulation  research  has  been  conducted 
by  Michigan  Tech,  MZA  Associates  Corporation,  ATK  Mission  Research,  Tempest  Technologies, 
and  Trex  Enterprises.  Independent  modeling  investigations  have  been  performed  by  these  member 
groups,  with  insights  from  all  influencing  the  wave-optics  propagation  models  developed  by  MZA 
Associates  Corporation  for  high-fidelity  simulation  of  directed-energy  weapons  systems.  The  MZA 
codes  have  been  used  extensively  by  UCLA,  Georgia  Tech,  and  Tempest  Technologies  for  investiga¬ 
tion  of  the  problems  of  adaptive  optics  and  target  tracking.  Thermal  blooming  models  developed 
by  Trex  Enterprises  have  been  incorporated  in  MZA’s  simulations  so  that  the  team  members  can 
develop  adaptive  optics  methods  and  image  processing  methods  to  mitigate  the  effects  of  thermal 
blooming  on  beam  control  and  tracking.  Also,  we  have  investigated  the  integration  of  the  novel 
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wavefront  sensing  methods  being  developed  at  Michigan  Tech  with  the  new  methods  for  adaptive 
optics  being  developed  at  UCLA  and  Tempest  Technologies. 

New  methods  for  correcting  wavefront  errors  have  been  developed  at  UCLA,  Tempest  Tech¬ 
nologies,  Michigan  Tech,  and  Trex.  The  effort  on  image  processing  methods  for  tracking  through 
turbulence,  the  primary  achievement  of  which  has  been  the  Bayesian  tracking  algorithm,  has  been  a 
collaboration  between  Georgia  Tech  and  Tempest  Technologies.  From  the  beginning  of  the  project, 
UCLA  and  Tempest  have  collaborated  closely  on  development  of  new  adaptive  optics  methods  based 
on  adaptive  filtering  and  control.  MZA  and  ATK  Mission  Research  have  worked  closely  with  UCLA, 
Georgia  Tech,  and  Tempest  to  provide  realistic  simulations  of  the  new  methods  for  beam  control  and 
tracking.  The  most  important  evaluations  of  the  new  adaptive  optics  methods  have  been  performed 
by  MZA  and  ATK  Mission  Research.  In  these  evaluations,  MZA  for  the  first  time  combined  one 
of  the  new  adaptive  optics  algorithms  developed  by  UCLA  with  the  Bayesian  tracker  developed 
by  Georgia  Tech  and  Tempest  Technologies.  That  these  new  methods  for  HEL  beam  control  and 
target  tracking  work  together  very  well  in  such  high-fidelity  simulations  is  one  of  the  greatest  success 
stories  of  this  MRI  to  date.  In  these  simulations,  the  adaptive  optics  algorithm  and  Bayesian  tracker 
worked  simultaneously  but  independently.  These  evaluations  by  MZA  suggest  that  planned  future 
research  to  integrate  the  two  methods  should  produce  even  greater  performance  enhancements. 

The  work  under  this  project  has  included  not  only  extensive  theoretical  and  numerical  research, 
but  also  substantial  experimental  research  at  Michigan  Tech,  UCLA,  MZA,  and  the  Air  Force  In¬ 
stitute  of  Technology  (AFIT).  UCLA  has  conducted  experimental  research  on  adaptive  control  of 
beam  jitter.  Also,  Michigan  Tech,  MZA,  ATK  Mission  Research  and  AFIT  have  conducted  ex¬ 
perimental  research  on  the  use  of  spatial  light  modulators  to  create  wavefront  errors  equivalent  to 
those  produced  by  atmospheric  turbulence.  Expanded  experimental  research  is  a  major  part  of  the 
research  planned  for  the  next  two  and  a  half  years  under  this  MRI.  The  experimental  research  will 
serve  two  purposes:  (1)  reveal  challenging  problems  inherent  in  the  physics  of  beam  control  and 
sensing,  thereby  guiding  our  theoretical  and  numerical  investigations  and  algorithm  development; 
(2)  demonstrate  the  practical  utility  of  the  methods  developed  in  our  research  for  beam  control  and 
tracking. 
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Chapter  2 

Adaptive  Control  of  Laser  Beams 

UCLA,  Tempest  Technologies, 

MZA  Associates  Corporation 


Compared  to  the  classical  methods  used  in  adaptive  optics  and  beam  steering  and  pointing,  the 
methods  developed  at  UCLA  can  improve  the  performance  of  directed-energy  weapons  such  as 
the  airborne  laser  (ABL)  and  the  airborne  tactical  laser  (ATL),  as  well  as  laser  communications 
systems.  High-fidelity  wave-optics  simulations  of  directed  energy  systems  have  shown  significant 
performance  improvements  with  the  new  adaptive  control  scheme  [l]1.  Also,  recent  experiments 
in  the  Atmospheric  Simulation  and  Adaptive-optics  Laboratory  Testbed  (AS ALT)  at  the  Starfire 
Optical  Range  at  the  Air  Force  Research  Laboratory,  Kirtland  AFB  have  demonstrated  enhanced 
performance  produced  by  UCLA’s  new  adaptive  control  scheme  for  adaptive  optics  [2,  3]. 

UCLA’s  contributions  to  control  of  laser  beams  fall  into  two  main  categories:  adaptive  control 
and  filtering  for  correction  of  higher-order  wavefront  errors  adaptive  optics  [1,  2,  3,  4,  5],  and  adaptive 
control  of  tilt  jitter  [6,  7,  8,  9,  10,  11,  12,  13,  14,  15,  16,  17,  18,  19].  The  two  problem  categories  are 
closely  related  and  are  both  present  in  adaptive  optics  systems.  Adaptive  control  of  higher-order 
wavefront  errors  usually  involves  control  loops  with  many,  often  hundreds  of  channels  corresponding 
to  spatially  distributed  wavefront  sensor  measurements  and  deformable  mirror  actuators;  adaptive 
control  of  tilt  jitter  involves  two  control  channels  and  two  or  more  sensor  signals,  but  usually  much 
higher  temporal  orders  of  the  adaptive  filter. 

A  brief  discussion  of  adaptive  optics  in  laser  weapon  systems  should  give  some  perspective  to 
much  of  the  proposed  research.  Adaptive  optics  (AO)  refers  to  the  use  of  deformable  mirrors  driven 
by  active  control  loops  that  feedback  wavefront  sensor  (WFS)  measurements  to  compensate  for 
turbulence- induced  phase  distortion  of  optical  waves  propagating  through  the  atmosphere.  These 
control  loops  reconstruct  (i.e.,  estimate  and  predict)  the  phase  profile,  or  wavefront,  from  the  WFS 
data.  The  control  loops  in  classical  AO  systems  are  linear  and  time-invariant  (LTI),  having  fixed 
gains  based  on  assumed  statistics  of  atmospheric  turbulence.  Such  control  loops  are  not  themselves 
adaptive  in  the  sense  in  which  the  term  adaptive  is  used  in  the  control  and  filtering  literature,  where 
adaptive  normally  refers  to  updating  control  and/or  filter  gains  in  real  time. 

Adaptive  compensation  is  needed  in  many  AO  applications  because  wind  velocities  and  the 
strength  of  atmospheric  turbulence  can  change  rapidly,  rendering  any  fixed-gain  reconstruction  al¬ 
gorithm  far  from  optimal.  UCLA  research  has  introduced  adaptive  wavefront  reconstruction  algo¬ 
rithms  that  use  recursive  least-squares  (RLS)  lattice  filters  to  predict  the  wavefront  and  estimate 
optimal  reconstructor  matrices  that  track  unknown  and  time-varying  turbulence  statistics.  In  this 

1UCLA  publications  can  be  found  at  http://www.beamcontrol.seas.ucla.edu. 
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approach,  an  adaptive  control  loop  augments  a  classical  AO  feedback  loop.  Results  in  [1,  20]  have 
shown  that  UCLA’s  adaptive  control  loops  are  robust  with  respect  to  modeling  errors  and  sensor 
noise. 

Figure  2.1  shows  a  schematic  diagram  for  an  adaptive  optics  problem  in  laser  weapon  system. 
This  diagram  represents  roughly  the  high-fidelity  ABL  simulation  in  [1] ,  as  well  as  an  array  of  laser 
weapons  motivating  the  AFRL  experimental  testbed  described  in  [2].  In  such  systems,  actuators 
are  distributed  in  a  two-dimensional  array  over  a  deformable  mirror.  These  actuators  are  driven  to 
adjust  the  profile  of  the  mirror  surface  and  cancel  the  phase  distortions  induced  in  the  high-energy 
laser  beam  as  it  propagates  through  atmospheric  turbulence.  A  wave  front  sensor  (WFS)  measures 
the  residual  wavefront  error,  using  an  array  of  subapertures  that  sense  the  spatial  derivatives,  or 
slopes,  of  the  phase  profile  on  a  grid  interlaced  with  the  locations  of  the  actuators. 

The  purpose  of  AO  system  is  to  compensate  the  outgoing  high  energy  laser  for  the  wavefront  error 
that  will  be  induced  by  atmospheric  turbulence,  so  that  the  laser  forms  a  fixed,  tight  spot  (image) 
on  the  target.  The  control  system  uses  a  beacon  created  by  illuminating  the  target  with  a  low 
energy  laser  as  the  basis  for  determining  the  commands  to  the  deformable  mirror  required  to  cancel 
turbulence-induced  phase  distortion.  Because  the  beacon  is  considered  to  be  a  distant  point  source, 
the  wavefront  propagating  from  the  beacon  would  be  very  nearly  a  plane  wave  when  it  reached 
the  mirror  with  no  atmospheric  turbulence.  This  plane  wave  is  the  desired  set  point  for  the  control 
algorithm.  If  the  wavefronts  propagating  from  the  beacon  to  the  target  travel  through  approximately 
the  same  atmosphere,  then  correcting  the  wavefront  from  the  beacon  should  compensate  for  the 
turbulence  effects  on  outgoing  beam. 


Figure  2.1:  Diagram  of  a  directed  energy  system  with  high  energy  laser  (HEL).  Key  components 
of  the  beam  control  system:  adaptive  optics  algorithm  (AO),  deformable  mirror  (DM),  wavefront 
sensor  (WFS). 


2.1  Collaboration  with  AFRL  and  Industry  on  Adaptive  Op¬ 
tics 

2.1.1  Recent  Results  and  Transitions  in  the  Testbed  at  the  Starfire  Op¬ 
tical  Range 

An  essential  feature  of  UCLA’s  research  supported  by  AFOSR  and  HEL  JTO  has  been  close  collabo¬ 
ration  with  the  Air  Force  Research  Laboratory  (AFRL)  and  industry.  UCLA  and  AFRL  researchers 
have  demonstrated  a  new  adaptive  control  scheme  for  adaptive  optics  in  experiments  in  the  Atmo- 
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spheric  Simulation  and  Adaptive-optics  Laboratory  Testbed  (AS ALT)  at  the  Starfire  Optical  Range 
at  the  Air  Force  Research  Laboratory,  Directed  Energy  Directorate,  Kirtland  AFB.  The  adaptive 
control  scheme  was  developed  at  UCLA  in  collaboration  with  researchers  at  the  Air  Force  Research 
Laboratory.  Initial  results  from  this  collaboration  were  presented  in  [2,  3]. 

Figure  2.2  shows  a  schematic  diagram  of  the  ASALT  optical  system.  As  shown  in  Figure  2.3,  the 
adaptive  control  loop  augments  a  classical  AO  loop  to  enhance  beam  control  and  imaging  through 
turbulence.  High-fidelity  wave-optics  simulations  of  directed  energy  systems  have  shown  signifi¬ 
cant  improvement  in  Strehl  ratio  (i.e.,  on-target  intensity)  and  tracking  jitter,  and  such  enhanced 
performance  now  has  been  confirmed  by  the  first  experimental  application  of  the  new  methods  [2] . 

Experiments  were  performed  in  the  ASALT  laboratory  to  evaluate  the  performance  of  the  quasi- 
adaptive  version  of  the  adaptive  control  loop.  In  the  experiments,  the  first  150  modes  from  a  set  of 
frequency-weighted  deformable-mirror  modes  were  used  by  the  adaptive  control  loop  [2].  First,  3000 
wavefront  sensor  frames  were  used  to  identify  the  adaptive  filter  gains,  and  then  the  performance  of 
the  adaptive  controller  was  evaluated  on  1000  frames  independent  of  those  used  for  identification. 

For  comparison,  the  same  experiment  was  performed  with  only  the  classical  AO  and  track  loops, 
using  the  same  1000  frames  for  evaluation.  For  the  turbulence  scenario  examined,  the  adaptive 
controller  provided  a  nearly  50%  increase  in  Strehl  ratio  and  reduced  the  variability  by  more  than 
15%.  Figure  2.4  shows  example  images  from  the  evaluation  sequences,  further  demonstrating  the 
benefits  of  the  adaptive  controller. 


2.1.2  Recent  Transition  to  Aero- Optics 

UCLA’s  methods  for  adaptive  control  in  adaptive  optics  are  being  used  in  a  Phase  II  SBIR  to  MZA 
Associates  Corporation  for  mitigation  of  aero-optics  effects  in  directed  energy  weapons,  funded  by 
MDA.  This  work  is  lead  by  Dr.  Matthew  Whitely  and  colleagues  at  MZA  Associates  Corporation, 
Albuquerque,  NM,  and  Dayton,  OH. 


2.1.3  Future  Collaboration 

Currently,  UCLA  researchers  are  working  with  AFR.L  researchers,  especially  Dr.  Darryl  Sanchez 
and  Lt.  Robert  Vincent,  at  the  Starfire  Optical  Range  (SOR),  Kirtland  AFB  to  implement  the 
fully  adaptive  version  of  UCLA’s  latest  adaptive  optics  algorithm  in  SOR’s  Atmospheric  Simulation 
and  Adaptive-optics  Laboratory  Testbed  (ASALT).  This  collaboration  should  produce  performance 
enhancements  in  the  ASALT  beam  control  system  and  should  lead  to  several  advances  in  adaptive 
optics.  Two  important  areas  of  research  will  be  the  design  of  adaptive  and  optimal  controllers  to 
mitigate  the  effects  of  wavefront  sensor  noise  in  adaptive  optics,  and  the  development  of  system 
identification  methods  for  obtaining  high-fidelity  models  of  complex  adaptive  optics  systems.  Both 
of  these  research  topics  will  address  critical  problems  that  significantly  limit  performance  in  beam 
control  systems  at  AFRL  and  in  most  defense  and  commercial  applications.  These  problems  can 
be  studied  to  some  extent  with  simulation  models  and  university  experiments;  however,  the  more 
realistic  testbed  at  AFRL’s  Starfire  Optical  Range  enables  much  more  relevant  research,  so  UCLA 
plans  to  continue  the  collaboration  represented  by  [2,  3] . 

UCLA  plans  to  collaborate  with  Teledyne  on  control  of  their  new  liquid  crystal  spatial  light 
modulator  for  wavefront  control  of  high  energy  lasers.  In  recent  months,  Professors  Gibson  and 
Tsao  have  been  invited  to  present  our  research  to  beam  control  groups  at  Northrop  Grumman  and 
Aerospace  Corporation.  Those  visits  initiated  plans  for  future  collaborations. 
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2.2  Recent  Results  and  Transitions  for  Control  of  Liquid 
Crystal  Beam  Steering  Devices 

UCLA,  AFRL  and  Teledyne  Scientific  Co.  have  collaborated  to  apply  feedback  and  adaptive  feedfor¬ 
ward  control  to  Teledyne’s  new  liquid  crystal  beam  steering  devices  [16,  15,  17].  These  novel  beam 
steering  devices  are  being  developed  as  actuators  in  jitter  control  and  adaptive  optics  for  applications 
to  laser  weapons  and  laser  communications.  Compared  to  standard  mirrors  used  for  beam  control, 
the  liquid  crystal  devices  have  the  advantages  of  low  power  consumption  and  no  moving  parts. 

Figure  2.5  and  2.6  show  UCLA’s  beam  control  experiment  with  the  prototype  liquid  crystal 
device.  Figure  2.7  shows  the  block  diagram  of  the  overall  control  system  including  the  liquid  crystal 
device,  the  optical  sensor,  the  disturbances  and  the  control  loops.  An  important  question  about 
the  new  devices  is,  can  they  deliver  the  control  bandwidths  comparable  to  those  of  fast  steering 
mirrors?  The  experimental  results  in  Figure  2.8  from  [17]  for  a  prototype  device  and  more  extensive 
experimental  results  in  [16,  15,  17]  are  quite  positive. 

UCLA’s  adaptive  control  loop  had  to  be  modified  to  handle  nonlinearities  in  the  liquid  crystal 
device  resulting  from  a  rate  limit  and  quantization  effects.  Such  nonlinearities  have  not  been  en¬ 
countered  in  our  research  with  fast  steering  mirrors  [21,  13,  22,  14,  23,  24],  but  our  newest  adaptive 
control  design  accommodates  the  nonlinearities,  resulting  in  no  apparent  performance  degradation. 

The  close  collaboration  among  UCLA,  AFRL  and  Teledyne  Scientific  Co.  has  been  very  produc¬ 
tive  in  several  ways.  First,  the  collaboration  has  illustrated  the  benefits  of  integrating  considerations 
of  control  system  performance  into  hardware  design.  Teledyne  based  the  re-design  of  the  driver  for 
the  two-axis  device  partly  on  the  performance  of  an  earlier  single-axis  device  in  control  experiments 
at  UCLA.  The  resulting  two-axis  device  allowed  the  adaptive  control  loop  to  achieve  much  higher 
bandwidths  than  with  the  initial  device.  Second,  UCLA  students  and  faculty  have  had  the  oppor¬ 
tunity  of  working  with  an  exciting  new  class  of  hardware  being  developed  in  industry  for  Air  Force 
missions.  The  experimental  results  reported  in  [16,  15,  17]  were  obtained  from  a  jitter  control  exper¬ 
iment  in  UCLA’s  beam  control  laboratory  with  the  Teledyne  liquid  crystal  device.  Most  recently, 
UCLA  Ph.  D.  student  Pawel  Orzechowski  has  worked  with  AFRL  and  Teledyne  researchers  to  set 
up  a  similar  experiment  at  the  Starfire  Optical  Range,  and  we  plan  to  continue  this  collaboration. 

There  are  three  nonlinearities  in  the  liquid  crystal  device:  a  rate  limit,  an  angle  saturation  and 
quantization.  The  effect  of  these  nonlinearities  is  investigated  in  detail  in  [17],  and  a  modification 
of  the  adaptive  controller,  based  on  a  nonlinear  model  of  the  plant,  is  introduced  to  compensate 
partially  for  the  nonlinearities.  The  analysis  in  [17]  shows  that  the  nonlinearities  are  significant  for 
the  jitter  levels  in  Figure  2.8,  but  the  original  UCLA  adaptive  controller,  which  does  not  take  the 
nonlinearities  into  account,  handles  the  effects  of  the  nonlinearities  as  well  as  the  modified  adaptive 
controller  does.  However,  as  illustrated  by  the  output  error  for  Axis  2  in  Figure  2.9,  when  the  jitter 
amplitudes  are  increased  by  only  33%,  the  adaptive  controller  based  on  a  linear  plant  model  fails, 
whereas  the  adaptive  controller  with  the  nonlinear  modification  performs  very  well.  These  results 
suggest  one  of  the  topics  of  the  proposed  research:  modeling  and  identification  of  nonlinearities 
in  liquid  crystal  beam  control  devices,  and  designing  controllers  that  take  the  nonlinearities  into 
account  to  optimize  the  performance  of  such  devices. 

Transition  to  Jitter  Control  in  Relay  Optics 

During  the  following  year,  UCLA’s  adaptive  jitter  control  methods  will  be  used  in  a  relay-optics 
experiment  at  AFRL  under  a  Phase  II  SBIR  to  Tempest  Technologies,  funded  by  MDA.  This  work 
is  lead  by  Dr.  Ben  G.  Fitzpatrick  of  Tempest  Technologies,  Los  Angeles,  CA. 
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Figure  2.2:  Atmospheric  Simulation  and  Adaptive-optics  Laboratory  Testbed  (ASALT)  optical  sys¬ 
tem,  Starfire  Optical  Range,  Kirtlancl  AFB.  The  Self-Referencing  Interferometer  Wavefront  Sensor 
(SRI  WFS),  an  innovative  sensor  being  developed  at  the  Starfire  Optical  Range,  was  used  for  the 
experiments  described. 


Figure  2.3:  Block  diagram  of  the  digital  control  loops  for  adaptive  optics.  ASALT  optics  block 
represents  the  optical  system  in  Figure  2.2. 


Classical  AO,  all  modes 


Adaptive  Control 


Figure  2.4:  Representative  closed-loop  scoring  camera  images.  Adaptive  control  produces  tighter 
laser  spot  with  greater  intensity  on  target  than  classical  adaptive  optics  (AO). 


Figure  2.5:  Left:  UCLA  beam  control  experiment  with  Teledyne’s  prototype  liquid  crystal  beam 
steering  device.  Control  sample  and  hold  rate  =  3125  Hz.  Right:  Close-up  view  of  the  liquid  crystal 
device. 


Laser  Source  Polarizer  Liquid  Crystal  Fast  Steering  Mirror  (FSM) 


Figure  2.6:  Diagram  of  UCLA  beam  control  experiment  with  Teledyne’s  liquid  crystal  beam  steering 
device. 


Figure  2.7:  Block  diagram  of  the  control  system,  ds  =  disturbance  command  to  shaker;  dg  = 
building  vibration;  dc  =  disturbance  command  FSM;  =  response  of  FSM;  0  =  beam  angle 
from  liquid  crystal  beam  steering  device;  y  =  beam  position;  OPS  =  optical  position  sensor;  u  and 
v  =  control  commands. 
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Figure  2.8:  Disturbance  rejection  performance  comparison  for  the  horizontal  axis  (Axis  1)  and  the 
vertical  axis  (Axis  2).  LTI  feedback  control  (red);  adaptive  control  (blue).  Maximum  lattice  filter 
order  N  =  10.  Jitter  sources:  building  vibration,  commands  to  fast  steering  mirror  (FSM)  and 
shaker. 
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Figure  2.9:  Performance  comparison  for  jitter-command  magnitudes  increased  by  33%:  LTI  feedback 
only,  adaptive  control  loop  with  linear  plant  model  G  and  adaptive  control  loop  with  nonlinear  plant 
model  Gml ■  Output-error  time  series  and  corresponding  power  spectral  densities  (PSD)  for  steady- 
state  responses  are  shown.  Axis  1  (left)  is  horizontal,  and  Axis  2  (right)  is  vertical.  Adaptive  control 
begins  at  t  =  19  sec. 
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Future  Collaboration 


UCLA  faculty  and  students  will  continue  the  current  collaboration  with  Dr.  Dan  Herrick  and  others 
at  AFRL  and  Teledyne  Scientific  Co.  on  control  of  liquid  crystal  devices  for  beam  steering.  While  the 
recent  experimental  research  in  UCLA’s  beam  control  laboratory  on  these  devices  has  been  quite 
productive,  it  is  even  more  exciting  now  that  UCLA’s  control  algorithms  are  being  implemented 
in  the  jitter  control  laboratory  at  AFRL’s  Starfire  Optical  Range  at  Kirtland  AFB.  Liquid  crystal 
technology  for  beam  control  also  is  being  developed  by  other  defense  contractors,  including  Raytheon, 
and  we  expect  to  develop  collaborations  with  such  companies. 
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Chapter  3 

HEL  Propagation  through 
Extended  Tubulence 

Michigan  Tech 


Michigan  Tech  has  conducted  significant  modeling  efforts  for  laser  beam  propagation  through  hori¬ 
zontal  paths  on  the  order  of  10  km  to  20  km  in  length  in  support  of  our  work  in  developing  advanced, 
higher  order,  nonlinear  beam  control  algorithms  for  deformable  mirrors.  This  work  has  been  aimed 
at  understanding  the  “forward  problem  of  beam  control  for  high  energy  lasers  in  a  tactical  scenario. 
Tactical  scenarios  differ  from  ABL-like  scenarios  in  that  the  propagation  path  is  at  lower  altitudes 
due  to  the  look-down,  shoot-down  aspect  of  tactical  engagements.  The  consequence  for  optical 
propagation  is  stronger  C\  levels  than  at  higher  altidues,  and  the  possibility  that  thermal  blooming 
effects  will  strongly  impact  the  system  performance.  The  isoplanatic  angle  is  expected  to  be  much 
smaller  than  the  field  of  view  of  the  target  acquisition  system  in  most  tactical  scenarios,  and  scin¬ 
tillation  is  expected  to  generally  be  strong.  Additionally,  there  will  not  be  a  point-like  beacon  in 
the  target  plane  for  tracking  and  higher  order  wave  front  sensing  purposes,  and  hence  one  must  be 
created  in  an  ABL-like  manner,  or  the  wave  front  information  must  be  extracted  from  the  scene. 
The  hardware  implications  for  tactical  HEL  systems  include  use  of  smaller  apertures  than  ABL-type 
systems,  and  the  use  of  lower  energy,  solid  state  lasers. 

The  physical  model  for  a  beacon  creation  system  in  a  tactical  HEL  is  shown  in  Fig.  3.1.  An 
uncompensated  outgoing  laser,  referred  to  as  the  beacon  laser,  is  used  to  illuminate  the  target. 
Since  there  is  no  signal  available  suitable  for  wave  front  sensing,  the  beacon  laser  is,  of  necessity, 
uncompensated.  As  a  result,  the  beam  arriving  at  the  target  is  affected  by  the  turbulence  present 
between  the  transmitter  and  the  target.  If  the  combination  of  path  length  and  turbulence  strength 
are  sufficiently  strong,  the  beam  arriving  at  the  target  will  be  on  average  broader  than  the  limit 
imposed  by  diffraction,  will  wander  randomly,  and  will  be  speckled  [25].  The  surface  of  the  target  is 
modeled  as  being  optically  rough  in  the  sense  that  the  fine  structure  of  the  scattering  surface  causes 
the  phase  of  the  scattering  surface  to  be  uniformly  distributed  on  (— 7r,7r)  [26].  Additionally,  the 
scattering  surface  is  (5-correlated  in  space,  and  its  rate  of  change  is  much  faster  than  the  integration 
time  of  wave  front  and  imaging  sensors.  At  the  telescope  pupil,  this  characteristic  of  the  target 
surface  will  give  rise  to  a  phenomenon  referred  to  as  laser  speckle  [26].  Laser  speckle  effects  include 
strong  intensity  and  phase  fluctuations  arising  from  the  random  phase  distribution  of  the  scattering 
surface.  Thus,  the  beacon  field  arises  from  a  random  intensity  incoherent  object  which  may  be  large 
relative  to  the  diffraction  limited  spot,  and  additionally,  at  the  wave  front  sensor  receiver  plane  the 
incoming  field  will  in  general  be  corrupted  with  laser  speckle  effects.  Some  of  the  radiation  scattered 
from  the  target  will  propagate  back  in  the  direction  of  the  transmit /receive  aperture,  where  it  will 
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Figure  3.1:  Physical  model  of  the  beacon  creation  process.  A  laser  passes  through  the  turbulent 
atmosphere  and  illuminates  the  target.  Backscattered  radiation  passes  back  through  the  atmosphere 
to  be  used  for  wave  front  sensing. 


be  intercepted  and  used  for  wave  front  sensing. 

The  beacon  field  can  be  measured  in  various  ways,  and  these  measurements  can  be  processed  to 
estimate  the  phase  of  the  impulse  response  for  propagation  through  the  atmosphere  from  the  aim 
point  on  the  target  to  the  laser  output  aperture.  This  impulse  response  is  denoted  by  Ha{xt,  xp,t), 
where  xp  is  a  coordinate  in  the  target  plane,  xp  is  accoordinate  in  the  laser  output  aperture  plane, 
and  t  represents  time.  In  general  hA(xT,  xp,t)  is  complex  valued: 


hA{xT,  xp,  t)  =  \hA{xT,  xp,  t)\  exp  \j(f>A(xT,  xP,  t)] 


(3.1) 


where  4>a(xt,  xp,  t)  is  the  phase  of  the  atmospheric  propagation  impulse  response,  and  under  aniso- 
planatic  conditions  Ha(xt ,  xp ,  t)  is  a  function  of  both  beacon  plane  position  and  pupil  plane  position. 
In  the  context  of  adaptive  optics  control  of  the  outgoing  laser  beam,  the  purpose  of  wave  front  sensing 
is  to  estimate  <Pa{xt,  Xp,t)  from  the  available  measurements.  It  should  be  noted  that  the  temporal 
behavior  of  turbulence  and  laser  speckle  effects  are  significantly  different  -  in  general  we  expect  that 
the  laser  speckle  effects  evolve  at  a  much  higher  rate  than  the  turbulence  effects,  and  this  insight 
must  be  incorporated  into  any  analysis  of  wave  front  sensing  in  this  environment. 

We  now  examine  the  characteristics  of  the  beacon  beam  arriving  at  the  target  plane.  We  shall 
begin  our  analysis  using  results  obtained  from  the  Rytov  theory  of  wave  propagation  through  tur¬ 
bulence,  which  assumes  weak  field  perturbations,  though  it  must  be  noted  that  the  assumption  of 
weak  perturbations  is  often  violated  in  practice  due  to  the  combination  of  turbulence  strength  and 
path  length.  We  shall  use  simulations  later  in  this  paper  to  obtain  results  for  conditions  where 
the  weak  fluctuation  assumption  is  violated,  with  the  result  that  the  variance  of  the  log-amplitude 
fluctuations  saturates  [25].  Consider  a  region  of  constant  C 2  over  a  path  of  length  L.  The  point 
statistics  of  the  intensity  fluctuations  of  the  beacon  beam  arriving  at  the  target  can  be  derived  from 
the  variance  of  the  log  amplitude  fluctuations  er^ .  For  a  collimated  beam  propagating  through  this 
a2  is  given  by  [25] 

a2  =  0.30 5As7/6C^i11/6,  (3.2) 

where  the  wave  number  is  given  by  =  27t/A,  where  A  is  the  wavelength.  The  normalized  variance 
of  the  intensity  fluctuations  at  the  target  erf  are  then  given  by 


a 


2 

i 


(V-Wf) 

W 


exp(4cr^)  -  1, 


(3.3) 


where  (•)  represents  the  statistical  expectation  opertor.  The  transverse  correlation  length  p0  of  a 


14 


collimated  beam  passing  through  this  turbulent  path  is  given  by  [25] 


p0  =  (lA6k2C2L)  3/5 . 

The  mean  square  short-term  beam  radius  due  to  turbulence  (p|)  is  [25] 


(P2s) 


4  L2 
JkDf 


fD\2  4  L2 

1-0.62  | 

tp0y/3' 

V  2  )  +  (kpo)2 

kd) 

-,6/5 


(3.4) 


(3.5) 


and  we  note  that  this  expression  describes  the  mean  square  instantaneous  spot  radius  in  the  target 
plane.  The  intensity  pattern  of  the  beacon  laser  at  the  target  will  also  wander  due  to  turbulence- 
induced  tilt,  leading  to  a  generally  much  larger  long-term  mean  square  spot  radius.  Since  we  are 
primarily  interested  in  making  high  speed  measurements,  evaluating  Eq.  (3.5)  is  sufficient  for  our 
purposes.  The  isoplanatic  angle  90  for  this  path  is  given  by  [25] 

90  =  (l.09 /c2C2L8/3)  3/5  ,  (3.6) 


and  the  isoplanatic  angle  projected  into  the  target  has  dimension  9qL. 

We  now  evaluate  ct2  and  both  ps  and  \9qL,  the  radius  of  a  spot  in  the  target  plane  subtending 
9q,  for  the  case  of  a  fixed  tactical  engagement-like  scenario,  with  path  length  L  =  20  km,  A  =  1  /jin, 
transmitter  diameter  of  D  =  1  m,  and  10-1'  <  C2  <  10~15  m-2/3.  We  have  two  purposes  for 
these  calculations:  (1)  to  discover  the  approximate  value  of  C2  for  which  a2  reaches  the  saturation 
regime  of  cr2  >  0.3  for  this  optical  path;  and  (2)  to  compare  the  short  term  spot  radius  to  the 
radius  of  a  0o~sized  spot  in  the  target  plane.  The  results  are  presented  in  Fig.  3.2.  Figure  3.2(a) 
presents  the  results  for  cr2  vs.  C2,  and  Fig.  3.2(b)  presents  the  results  for  ps  and  \90L  at  the  target 
plane  vs.  C2.  Inspection  of  Fig.  3.2(a)  shows  that  for  the  atmospheric  path  modeled,  turbulence 
characterized  by  C2  <  2  x  10-16  yields  values  of  cr2  which  are  below  the  saturation  regime.  At 
altitudes  below  2  km  C2  values  in  the  range  10-17  <  C2  <  10-15  m-2/3  are  common,  and  we 
conclude  that  many  scenarios  of  practical  interest  for  low  altitude  beam  projection  systems  will 
result  in  the  need  to  operate  in  the  presence  of  saturated  cr2  conditions  [25] .  Higher  altitude  systems 
will  also  experience  a  similar  operational  situation,  though  longer  optical  paths  are  permitted  before 
saturated  cr2  conditions  arise. 

Inspection  of  Fig.  3.2(b)  shows  that  the  Rytov  theory  predicts  that  the  spot  radius  due  to  short 
term  beam  spreading  is  in  the  range  of  1.8  to  more  than  20  times  the  radius  of  an  isoplanatic  angle¬ 
sized  patch  in  the  target  plane,  \90L.  If  we  assume  the  target  is  on  the  order  of  the  same  size  as 
ps  plus  the  beam  wander,  or  bigger,  then  most  or  all  of  the  light  arriving  at  the  target  plane  will 
be  scattered.  If  the  target  is  optically  rough,  some  light  from  each  illuminated  point  on  the  target 
plane  will  be  scattered  back  in  the  direction  of  the  aperture.  Because  ps  \9qL,  light  arriving 
at  the  aperture  from  different  points  on  the  target  will  have  travelled  through  significantly  different 
atmospheric  paths.  Light  corrupted  with  the  aberrations  obtained  from  these  many  optical  paths 
arrives  superimposed  at  the  aperture,  where  it  is  used  in  the  wave  front  sensing  system.  Because 
the  light  arrives  superimposed  from  so  many  directions  which  are  in  general  separated  by  more  than 
f?o,  conventional  wave  front  sensing  processes  will  not  estimate  the  same  turbulence-induced  wave 
front  error  that  would  be  computed  if  a  point  source  beacon  were  present.  In  most  cases  of  practical 
interest,  the  errors  between  the  wave  front  estimated  from  the  extended  beacon,  and  the  wave  front 
which  would  have  been  estimated  if  a  point  source  beacon  were  present  are  significant,  and  severely 
degrade  the  ability  of  the  beam  control  system  to  focus  light  on  a  small  spot  in  the  target  plane 
[27].  This  is  the  essence  of  the  effect  referred  to  as  beacon  anisoplanatism.  This  topic  is  addressed  in 
more  detail  in  the  paper  Fundamental  considerations  for  wave  front  sensing  with  extended  random 
beacons,  by  Roggemann,  which  was  presented  at  the  SPIE  meeing  in  Denver,  in  August,  2004.  This 
paper  is  included  in  the  appendix. 
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(a)  (b) 

Figure  3.2:  (a)  a £  at  the  target  plane  vs.  C\  for  L  =  20  km,  and  A  =  1  /rm.  (b)  ps  and  \9qL  at 
the  target  plane  vs.  for  the  same  optical  path,  with  transmitting  aperture  diameter  fl  =  lm. 


The  key  findings  of  the  modeling  work  described  above  are  that  the  total  implications  of  the 
lack  of  a  cooperative  beacon  in  tactical  HEL  scenarios  is  likely  the  most  fundamental  performance 
limitation.  Both  tracking  and  higher  order  compensation  are  affected  by  this  issue.  As  we  will  show 
later,  when  an  appropriate  tracking  signal  is  present,  advanced,  non-linear  deformable  mirror  control 
algorithms  are  very  effective  in  providing  a  high  Strehl  ratio  at  the  target.  Likely,  technological 
solutions  to  the  these  problems  for  HEL  beam  control  will  have  the  following  key  elements: 

1.  Tracking  information  obtained  from  some  element  of  the  target  or  scene.  Work  along  the  lines 
of  the  extended  target  tracking  efforts  underway  at  Georgia  Tech  offers  hope  in  this  area,  but 
the  algorithms  will  have  to  be  extended  to  work  in  a  cluttered  background  environment. 

2.  A  beacon  laser  will  still  be  required  to  provide  higher  order  wave  front  sensing  information,  and 
processing  the  resulting  measurements  into  useful  deformable  mirror  commands  will  require 
non-linear  algorithms. 

3.  Extremely  high  speed  computing  will  be  required  for  these  approaches  to  be  successful. 

We  now  move  on  to  a  top  level  discussion  of  the  extensive  simulation  efforts  underway  at  Michigan 
Tech. 

Simulation  of  a  tactical  HEL  path  requires  a  multi-layer  model  of  the  turbulence  to  capture 
the  scintillation  and  anisoplanatic  aspects  of  the  turbulence  effects.  Since  the  optical  path  we  have 
modeled  so  far  has  constant  C%  we  model  the  turbulence  with  10  equal  strength,  but  statistically 
independent  phase  screens.  These  screens  are  generated  with  a  von  Karman  power  spectral  density 
using  a  well-established  technique  based  on  spatial  filtering  of  a  white  noise  process.  The  first  screen 
is  placed  in  the  aperture  of  the  transmit/receive  telescope,  and  the  last  is  placed  a  finite  distance 
in  front  of  the  target.  Each  screen  is  treated  as  a  phase  object,  so  that  passage  through  the  screen 
causes  the  phase  of  the  optical  wave  to  change,  but  not  its  amplitude.  A  wave  optics  propagator 
based  on  the  angular  spectrum  propagator  [28]  is  used  to  model  all  of  the  propagations  between 
screens,  and  between  the  last  screen  to  the  target.  This  approach  has  been  widely  used  in  many 
propagation  studies,  and  is  the  one  implemented  in  widely  used  programs  such  as  WaveTrain.  We 
chose  the  approach  because  it  was  very  straight  forward  to  extend  the  vast  library  of  existing  code 
we  had  developed  and  validated  under  other  programs,  and  because  it  offered  maximum  flexibility 
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in  developing  the  non-linear  algorithms  for  deformable  mirror  control.  Typical  results  of  simulations 
of  this  sort  are  shown  elsewhere  in  this  report,  and  our  results  are  very  similar  to  these. 
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Chapter  4 

Novel  Methods  for  Wavefront 
Sensing  and  Beam  Control 

Michigan  Tech 


4.1  Optimal  Beam  Control 

We  consider  beam  propagation  scenarios  for  which  a  pupil-plane  field  f{u)  over  an  aperture  region 
u  €  A  is  transmitted  through  a  propagation  medium  that  is  characterized  by  an  inhomogenous 
and  random  medium  with  an  associated  Green’s  function  h(x,u),  where  a  is  a  two-dimensional 
spatial  index  in  the  pupil  plane  and  a:  is  a  two-dimensional  spatial  index  in  the  target  plane.  The 
target-plane  field,  then  is 

g{x)=  /  h(x,  u)f(u)du, 

Ja 

which  is,  in  general,  a  random  quantity  because  of  the  randomness  in  the  propagation  kernel  h.  The 
objective  of  wavefront  sensing  and  beam  control  is  to  generate  a  transmitted  beam  that  maximizes 
in  some  sense  the  localized  energy  |g(a;)|2  at  the  target.  One  common  method  for  quantifying  the 
localized  energy  of  a  transmitted  beam  is: 


Ir 


\g(x)\2dx, 


so  that  Ir  quantifies  the  total  intensity  over  a  disc  of  radius  r.  Because  the  propagation  kernel  is 
random,  the  localized  intensity  Ir  is,  in  general,  a  random  quantity. 


4.2  Time-averaged  optimal  beam  control 

One  approach  to  dealing  with  the  randomness  of  the  propagation  kernel  is  to  optimize  first-  or 
second-order  moments  of  the  localized  intensity  Ir.  We  have  shown  that  when  one  chooses  to 
optimize  the  expected  intensity  E[Ir],  then  the  optimal  beam  is  the  principle  eigenfunction  (the 
largest  eigenvalue’s  eigenfunction)  for  the  following  kernel: 

H(u,u')  =  I  E  [h* (x,u)h(x,u')\  dx. 

J  |rr|<r 
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Initial  studies  have  shown  that,  for  random  media  such  as  short-path  or  extended-path  turbulence, 
a  beam  specified  in  this  manner  provides  little,  if  any,  performance  enhancement  over  an  uncom¬ 
pensated  beam.  Because  of  this,  beam  control  must  be  accomplished  for  each  individual  realization 
of  the  random  Green’s  function. 


4.3  Optimal  beam  control  when  the  propagation  kernel  is 
known 

Optimization  of  the  localized  energy  of  a  transmitted  beam  through  a  medium  that  is  characterized 
by  the  Green’s  function  h{x ,  u)  requires  that  the  transmitted  beam  be  determined  as  the  principle 
eigenfunction  for  the  following  kernel; 

7i(u,u')  =  /  h*(x,u)h(x,u')dx. 

J \x\<r 

We  have  also  shown  that  this  eigenfunction  can  be  efficiently  computed  using  an  iterative  transform 
algorithm,  and  this  method  provides  garaunteed  convergence  to  the  optimal  beam.  However,  the 
Green’s  function  is  rarely  known,  so  some  method  of  adaptive  inference  must  be  utilized  in  most 
situations. 


4.4  Field  Models  for  Noncooperative  Targets 

Wavefront  sensing  and  beam  control  in  the  presence  of  a  noncooperative  target  typically  requires 
that  the  target  itself  be  illuminated  to  form  a  reference  field  for  sensing  and  control.  The  field  that 
presents  itself  at  the  sensor’s  pupil,  then,  can  be  used  to  estimate  the  wavefront  and  make  inferences 
about  the  propagation  medium.  However,  because  the  reference  field  that  reflects  from  the  target 
has  extended  spatial  extent,  the  field  at  the  sensor  pupil  does  not  allow  for  direct  inference  about 
the  medium. 

We  refer  to  the  complex  amplitude  of  the  field  that  is  reflected  from  the  target  as  gt{x),  and  this 
field  is  related  to  the  field  that  appears  in  the  pupil  as: 


fp(u)  =  j  h*(x,  u)gt{x)dx, 


where  we  have  assumed  that  reciprocity  holds  for  the  the  medium’s  propagation  kernel.  That  is,  the 
kernel  for  propagation  from  the  pupil  to  the  target  is  h(x,  u),  and  the  kernel  from  propagation  from 
the  target  to  pupil  is  h*(x,u).  Whereas  we  assume  that  the  medium  is,  in  general,  inhomogeneous 
and  its  Green’s  function  is  unknown,  we  do  assume  that  the  Green’s  function  does  not  change  during 
a  particular  observation  interval.  The  mutual  intensity  for  the  pupil-plane  field  over  an  observation 
interval  is: 

Jp{u,u’)  =<  fp(u)f*(u ')  >, 

where  <  •  >  denotes  a  time  average  over  the  observation  interval.  Accordingly, 

Jp(u,u')  =  <f(u)f*(u')> 

=  J  J  h*(x,u)  <  gt(x)gl(x')  >  h(x' ,  u')dxdx' 

=  J  J  h*(x,u)Jt(x,x')h(x',u')dxdx', 

where  Jt(x,x')  is  the  mutual  intensity  for  the  field  reflected  from  the  target. 
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uncompensated  beam 


compensated  beam 


Figure  4.1:  An  example  illustrating  pupil-plane  adaptive  compensation  with  a  noncooperative  target. 
Both  the  uncompensated  and  compensated  beam  profiles  are  shown. 

If  the  observation  time  interval  is  long  compared  with  the  coherent  speckle  variations  of  the 
reflected  field,  then  the  target-plane  mutual  intensity  will  be  well  approximated  by: 

Jt.(x,x')  ~  b(x)6(x,x'), 

where  b(x)  is  the  average  intensity  that  is  reflected  at  the  target,  sometimes  called  the  ’beacon’. 
Accordingly,  the  mutual  intensity  in  the  pupil  plane  is 

Jp(u,u/)~  J  h*(x,u)h(x,u')b(x)dx. 

4.5  Pupil-plane  Wavefront  Estimation  and  Beam  Control 

The  principle  eigenfunction  for  the  pupil-plane  mutual  intensity,  when  used  to  propagate  a  beam  to 
the  target,  will  have  the  property  that 

Ib  =  J  b(x)\g(x)\1 2 3dx 

is  maximized.  Because  of  this,  this  eigenfunction  can  be  used  in  an  adaptive  optics  system  for 
’optimal’  beam  control.  The  measurement  and  processing  steps  are  outlined  below: 

1.  Utilize  a  shearing  interferometer  or  similar  instrument  to  measure  the  pupil-plane  mutual 
intensity; 

2.  Determine  the  principle  eigenfunction  for  the  pupil-plane  mutual  intensity;  and 

3.  Utilize  the  principle  eigennfunction  as  the  transmitted  beam. 

An  example  of  the  utilization  of  this  method  is  shown  in  Figure  4.1.  The  beacon  intensity,  along 
with  the  uncompensated  and  compensated  beams  are  shown  in  Figure  4.2. 

The  difficulty  with  this  method  is  that  the  pupil-plane  mutual  intensity  is  a  four-dimensional 
function,  and  its  measurement  can  be  challenging.  In  the  next  section  we  discuss  a  method  whereby 
image-plane  measurements  are  used  to  estimate  and  correct  the  field. 
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Figure  4.2:  The  beacon  intensity  and  beam  intensities  for  the  pupil-plane  adaptive  compensation 
example  with  a  noncooperative  target. 


4.6  Image-plane  Wavefront  Estimation 


Suppose  an  adaptive  field-compensator  is  placed  in  the  system’s  pupil,  and  an  image-plane  intensity 
is  recorded  in  the  focus  plane  for  the  system’s  optics.  This  image  intensity  will  be  related  to  the 
pupil-plane  field  as: 

„  2 


i(y) 


fp{u)^{u)e 


-jllL 


u-y 


du 


where  i/j(u )  =  a(u)e is  the  pupil-plane  field  compensation.  If  the  compensation  is  selected  to 
maximize  the  image  sharpness  defined  as: 


5  = 


s(y)i(y)dy , 


where  s{y)  is  a  non-negative  sharpness  window,  then  the  optimal  compensation  will  be  the  principle 
eigenfunction  for  the  following  kernel: 

Hs(u,u')  =  J  h* (x,u)h(x,  ul)b(x)dxS  ^  ^ 

where  S  is  the  Fourier  transform  for  the  sharpness  window.  An  example  of  the  utilization  of  this 
method  with  a  Gaussian  sharpness  window  with  a  full-width  at  half  maximum  (FWHM)  equal  to 
Adi/D  is  shown  in  Figure  4.3.  The  beacon  intensity,  along  with  the  uncompensated  and  compensated 
beams  are  shown  in  Figure  4.4. 
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uncompensated  beam 


Figure  4.3:  An  example  illustrating  image-plane  adaptive  compensation  with  a  noncooperative  tar¬ 
get.  Both  the  uncompensated  and  compensated  beam  profiles  are  shown. 
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Figure  4.4:  The  beacon  intensity  and  beam  intensities  for  the  image-plane  adaptive  compensation 
example  with  a  noncooperative  target. 
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Figure  4.5:  Deformable  mirror  control  based  on  optimizing  an  image  sharpness  metric. 


4.7  Non-linear  Higher  Order  Deformable  Mirror  Control 

Significant  performance  limitations  arise  for  wave  front  sensing  and  deformable  mirror  control  in 
the  tactical  HEL  environment  from  beacon  anisoplanatism,  strong  scintillation,  and  the  lack  of  a 
cooperative  tracking  beacon.  Phase  difference  measuring  devices,  such  the  Hartmann  sensor,  are 
subject  to  large  errors  from  these  issues.  Having  experimented  extensively  with  phase  difference 
measuring  devices  in  the  tactical  HEL  scenario,  we  decided  to  pursue  a  deformable  mirror  control 
scheme  which  does  not  strictly  work  in  phase  space,  but  rather,  seeks  to  maximize  some  intensity- 
based  measure  of  performance  by  controling  the  phase  falling  on  the  pupil  with  a  deformable  mirror. 
Our  work  is  an  extension  of  an  idea  which  has  its  basis  in  the  image  sharpness  metric  work  of  Muller 
and  Buffington  [29],  and  has  been  recently,  and  successfully  extended  by  Vorontsov  [30].  Deformable 
mirror  control  algorithms  based  on  these  image  quality  measures  are  inherently  non-linear,  since 
there  is  no  linear  relationship  between  the  phase  correction  and  the  intensity  measurement.  We 
have  shown  under  this  project  that  deformable  mirror  control  based  on  this  concept  can  provide  a 
significant  performance  improvement  over  conventional  Hartmann  sensor-base  wave  front  control, 
but  at  the  cost  of  requiring  significant  real-time  processing  capability.  Excellent  tracking  of  the 
target  to  remove  turbulence-induced  tilt  errors  must  be  provided  in  either  case. 

The  block  diagram  for  the  system  we  are  studying  is  shown  in  Fig.  4.5.  An  uncompensated 
beacon  laser  is  used  to  illuminate  the  target.  Light  scattered  from  the  target  is  received  by  the 
aperture,  and  an  image  is  formed.  The  deformable  mirror  control  algorithm  seeks  a  set  of  deformable 
mirror  commands  which  optimize  an  image  sharpness  metric  in  the  image  plane.  When  the  optimal 
deformable  mirror  figure  is  found,  the  HEL  is  fired.  In  simulation  the  HEL  beam  is  propagated  to 
the  target,  and  various  measures  of  performance  are  computed.  It  must  be  noted  that  this  approach 
requires  tilt  control  information  be  obtained  from  some  source  other  than  the  beacon  laser,  since 
due  to  reciprocity,  the  return  beam  has  no  tilt  information.  For  our  purposes,  we  presume  that  a 
tilt-only  beacon  is  present  in  the  scene,  and  only  use  the  light  emanating  from  this  source  to  obtain 
tilt  correction  information.  In  an  actual  application  information  from  the  scene  itself  would  likely 
be  used  to  get  tilt  correction  information. 

The  basic  idea  of  this  deformable  mirror  concept  is  to  treat  the  operation  of  finding  the  optimal 
deformable  mirror  commands  as  an  optimization  problem.  Let  the  deformable  mirror  command 
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Figure  4.6:  (a)  Target  intensity  due  optimizing  the  integrated  square  intensity  Ji(a);(b)  associated 
target  intensity  which  would  arise  with  only  tilt  correction.  In  both  cases 
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vector  be  represented  by  a.  In  the  paradigm  describe  above,  the  measured  image  of  the  scattered 
beacon  laser  is  a  function  of  a  which  we  refer  to  as  I(xi,a),  where  Xi  is  an  image  plane  coordinate. 
Our  presumption,  which  has  been  born  out  by  simulation  experiments,  is  that  the  a  which  optimizes 
a  measure  of  the  quality  of  I(xi,  a)  is  an  appropriate  correction  to  apply  to  the  outgoing  beam.  We 
iterate  on  the  vector  a  to  optimize  an  objective  function  we  will  generically  refer  to  Jm(a),  where 
the  subscript  m  refers  to  an  index  on  a  specific  image  quality  metric.  As  an  example,  consider  use 
of  the  Ji(a)  metric  proposed  by  Muller  and  Buffington  given  by 

J\{a)  =  J  I2(xi,a)dxi  (4-1) 

Figure  4.6  shows  the  result  of  one  optimization  of  Ji(a)  for  a  tactial  HEL  scenario  with  C 2  = 
9.56  x  10-16m-2/3,  with  tilt  correction  present.  Figure  4.6(a)  shows  the  result  for  propagation 
through  a  single  realization  of  an  extended  atmosphere  in  tactical  HEL-like  conditions  using  adaptive 
optics  compensation;  Fig.  4.6(b)  shows  the  results  for  the  case  where  only  tilt  correction  was  present. 
Comparison  of  these  two  results  shows  clearly  that  deformable  mirror  control  based  on  non-linear 
optimization  of  image  sharpness  metric  holds  significant  promise. 

While  some  success  has  been  achieved  using  the  Ji(a)  metric,  we  have  obtained  better  perfor¬ 
mance  with  metric 

Jz{at)  =  f  I(xi,a)M(xi)dxi,  (4.2) 

where  M(xi)  is  a  mask  in  the  image  plane.  Essentially,  the  Js(a)  metric  seeks  to  find  a  deformable 
mirror  command  vector  a  which  maximizes  the  energy  falling  inside  some  region  in  the  image 
detector  plane.  We  have  had  excellent  results  with  this  approach,  and  an  example  is  presented  in 
Fig.  4.7.  Figure  4.7  shows  the  average  encircled  energy  in  the  target  plane  for  the  tactical  HEL 
scenario  with  C2  =  19.12  x  10-16m-2/3.  The  three  curves  shown  are  encircled  energy  for  the  cases 
of  no  compensation  (the  NOCOMP  curve),  tilt  correction  using  a  cooperative  beacon  plus  higher 
order  correction  based  on  least  squares  reconstruction  from  Hartmann  wave  front  sensor  outputs 
(the  HWFS  curve),  and  for  tilt  correction  using  a  cooperative  beacon  plus  higher  order  correction 
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Figure  4.7:  Percentage  encircled  energy  in  the  target  plane  for  C2  =  19.12  x  10  16m  2//3. 


obtained  from  optimizing  the  J3  metric  (the  NLOPT  curve).  Inspection  of  Fig.  4.7  The  promise  of 
the  non-linear  approach  to  controlling  the  deformable  mirror  is  clear.  The  complete  set  of  results 
obtained  so  far  are  included  in  the  Kizito  dissertation  provided  as  an  addendum  to  this  report. 
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Chapter  5 


Modeling  and  Simulation  of  HEL 
Wavefront  Propagation 

MZA  Associates  Corporation 

5.1  Introduction 

Since  April  of  2006,  MZA  Associates  Corporation  (MZA)  has  been  supporting  the  University  of 
California  at  Los  Angeles  (UCLA)  as  a  subcontractor  on  a  grant  from  the  Air  Force  Office  of 
Scientific  Research  (AFOSR)  sponsored  by  the  High  Energy  Laser  Joint  Technology  Office  (HEL- 
JTO).  This  report  summarizes  the  work  accomplished  over  the  period  from  1  April  2006  through 
31  October  2007.  MZA  supported  this  effort  in  the  following  areas: 

•  Phase  screen  generation  for  numerical  simulation  of  long  time  series  laser  beam  propagation 
through  atmospheric  turbulence. 

•  Closed  loop  tracking  performance  as  a  function  of  the  spatial  frequency  content  of  atmospheric 
turbulence  screens. 

•  Matlab®  version  of  the  North  American  Aerospace  Defense  Command  (NORAD)  SGP4/SDP4 
orbital  propagator  and  study  on  energy  transfer  to  orbital  assets. 

•  Implementation  of  the  Extreme  and  Percentile  Environmental  Reference  Tables  (ExPERT) 
atmospheric  data  into  the  Scaling  for  HEL  and  Relay  Engagements  (SHaRE)  Matlab®  toolbox. 

•  Implementation  of  an  adaptive  controller  (developed  by  UCLA)  in  Matlab®/ Simul  ink®  and 
WaveTrain™. 

•  Laboratory  simulation  of  atmospheric  turbulence  and  aero-optic  disturbances  using  spatial 
light  modulators. 

5.2  Accomplishments 

Over  the  19  months  of  this  effort,  MZA  has  accomplished  a  significant  amount  of  research  on  high 
energy  laser  (HEL)  propagation  through  the  atmosphere.  The  majority  of  the  work  was  related  to 
the  generation  and  implementation  in  WaveTrain™  of  phase  screens  for  long-time  series  closed-loop 
wave  optics  simulations.  Much  of  this  work  was  an  extension  of  the  work  completed  under  this  effort 
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by  Mission  Research  Corporation  and,  subsequently,  Alliant  Techsystems  (ATK).  A  brief  summary 
of  the  each  of  the  accomplishments  listed  above  will  be  given  in  this  section. 


5.2.1  Phase  Screens  for  Long  Time  Series  Wave  Optics  Simulations 

In  this  section  we  describe  the  progress  made  towards  implementation  of  new  atmospheric  phase 
screens  in  WaveTrain™  which  has  been  made  by  integrating  the  products  of  a  series  of  separate 
efforts  within  MZA.  The  overall  effort  has  been  subsidized  by  several  funding  sources  including  in¬ 
ternal  MZA  funding,  Small  Business  Innovative  Research  (SBIR)  funding,  the  Airborne  Laser  (ABL) 
program,  as  well  as  this  contract.  The  overall  effort  has  amounted  to  developing  new  atmospheric 
modeling  classes  within  the  WaveTrain™  code.  The  following  subsections  summarize  the  basic  im¬ 
provements  which  have  been  implemented  as  a  result  of  this  effort 

Atmospheric  Modeling  Classes 

When  the  atmospheric  model  for  WaveTrain™  was  originally  developed  in  1996,  it  was  implemented 
in  a  manner  so  as  to  achieve  certain  efficiencies  which,  although  appropriate  for  computational 
resources  at  the  time,  have  now  become  irrelevant  due  to  increased  processor  speed  and  memory. 
As  a  result,  the  code  was  not  as  modular  and  was  made  more  complex  than  it  now  needs  to  be. 
One  goal  of  this  effort  is  to  simplify  and  generalize  the  interface  to  the  atmospheric  model  to  make 
it  more  expandable  and  flexible  so  that  new  capabilities  can  be  easily  added  using  MZA’s  tempus™ 
plug-and-play  methodology. 

The  new  atmospheric  model  is  implemented  as  a  series  of  new  C++  classes  integrated  into 
WaveTrain™.  Descriptions  of  the  major  new  WaveTrain™  classes  are  provided  in  Table  5.1.  The 
wtPhaseScreen  class,  and  those  derived  from  it,  are  documented  with  inline  comments  in  the  Wave- 
Train™  source  code  repository.  The  PhaseScreen  and  Atmosphere  classes  are  also  documented  inline 
and  provide  standard  WaveTrain™  system  documentation  in  the  WaveTrain™  System  library. 


Table  5.1:  Descriptions  of  the  major  new  atmospheric  WaveTrain™  classes 


Class  Name 

Description 

wtPhaseScreen 

The  base  class  for  the  implementation  of  different  phase 
screen  computational  and  representational  models.  All 
new  phase  screen  implementations  derive  from  this  base 
class  to  allow  the  wrapper  code  contained  in  the  following 
two  classes  to  manipulate  different  phase  screen  implemen¬ 
tations  with  the  same  logic 

PhaseScreen 

tempus™  System  which  implements  a  phase  screen  logic 

Atmosphere 

Composite  tempus™  system  which  wraps  logic  around  a 
sequence  of  PhaseScreen  objects  in  order  to  implement 
the  complete  atmosphere 

New  Method  for  computing  Phase  Screens 

The  original  phase  screens  implemented  in  WaveTrain™,  based  on  a  Fourier  spectrum  method,  can  be 
made  to  approximate  correct  low  spatial  frequency  statistics  at  a  sometimes  significant  computational 
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cost.  MZA  has  implemented  two  new  methods  for  generating  phase  screens  which  are  approximately 
correct  at  low  spatial  frequencies. 

A  major  drawback  of  the  Fourier  spectrum  method  of  generating  phase  screens  is  that  it  requires 
the  phase  screen  to  be  generated  over  a  finite  region  predefined  at  the  time  the  screen  is  generated. 
This  results  in  having  to  store  in  memory  the  entire  phase  screen  all  at  once  even  when  only  a 
small  segment  is  needed  and  is  very  inefficient  in  the  cases  when  the  low  frequency  statistics  of 
the  screen  are  important.  The  two  new  phase  screen  generation  methods  are  potentially  infinitely 
extensible  in  order  to  allow  slewing  through  very  large  areas  of  the  phase  screen  without  having 
to  have  held  the  entire  screen  in  memory  all  at  once.  Furthermore,  this  capability  provides  the 
possibility  that  disconnected  segments  of  the  same  correlated  screen  can  be  generated  to  handle 
multi-static  propagation  cases. 

Previously,  Dr.  Eric  Magee  developed,  under  this  effort  and  while  at  Mission  Research  Corpora¬ 
tion  and  ATK,  a  technique  for  generating  infinitely  extensible  phase  screens  having  approximately 
correct  low  spatial  frequency  statistics  using  a  Fourier  series  (FS)  approach.  In  this  method,  a  series 
of  random  Fourier  coefficients  are  generated  whose  sum  when  evaluated  at  a  given  point  within  the 
phase  screen  results  in  an  optical  path  difference  (OPD)  having  the  correct  magnitude  and  rela¬ 
tionship  to  neighboring  points  to  represent  an  arbitrarily  specified  power  spectral  density  (PSD)  or 
structure  function.  For  atmospheric  phase  screens  the  PSD  used  is  typically  Kolmogorov,  however, 
the  same  routine  has  been  used  to  generate  Gaussian-correlated  phase  screens  having  similar  char¬ 
acteristics  to  aero-optical  effects.  This  capability  was  recoded  in  C++  and  was  added  to  WaveTrain™ 
as  the  class  FSScreen  which  is  derived  from  the  base  class  wtPhaseScreen. 

The  technique  is  based  on  the  decomposition  of  the  random  phase  into  a  modified  FS  expansion, 

N/2  N/2 

<t>{x,y)  =  E  E  Cn,m  exp  \j2n  (fXnx  +  fVmy)]  (5.1) 

n=—N/2m=  —  N/2 


where  the  are  the  random  complex  coefficients,  and  fxn  and  fym  are  the  user  defined  x—  and 
y— spatial  frequencies.  The  Frequencies  are  not  constrained  to  harmonics  as  in  the  discrete  Fourier 
transform  (DFT).  We  use  logarithmically  spaced  frequencies 
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where  the  fn  are  the  spatial  frequencies,  a  is  the  logarithmic  growth  parameter,  and  N  is  the  number 
of  frequencies  between  the  minimum  spatial  frequency,  /mjn,  and  the  maximum  spatial  frequency, 
/max-  The  statistics  of  the  cHjm  are  circular  complex  Gaussian  with  variance  determined  by  the 
desired  PSD  of  the  phase  screen 

E  {\cn,m\2}  =  AfXnAfVm<f>k(fXn,fyJ  (5.6) 

where  where  &k(fx,fy)  is  the  PSD  of  the  phase  variations.  The  main  advantage  of  the  FS  phase 
screen  is  that  once  the  coefficients  and  spatial  frequencies  are  determined,  the  phase  screen  is  defined 
for  all  points  in  a  plane.  An  example  phase  screen  at  various  resolutions  is  shown  in  Fig.  5.1(a)-(c). 
The  sample  phase  structure  function  is  shown  in  Fig.  5.1(d).  This  plot  also  shows  the  theoretical 
structure  function  (in  black)  and  the  expectation  of  the  structure  function  given  the  spatial  frequency 
content  of  the  phase  screen. 
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Figure  5.1:  Example  phase  screens  [(a)-(c)]  generated  using  a  modified  FS  with  a  minimum  spatial 
frequency  of  0.0001  m  .  The  structure  function  for  an  ensemble  of  5000  independent  screens  is 
shown  in  (d). 


Under  a  separate  effort,  Dr.  Russ  Butts  developed  a  technique  for  correct  and  extensible  phase 
screens  based  on  the  conditional  probability  density  functions  (CPDF).  The  original  implementation 
was  in  Matlab®.  Dr.  Richard  St.  John  developed  a  Java-based  tool  for  constructing  the  screens 
in  a  manner  consistent  with  the  wave-optics  codes.  As  shown  in  Fig.  5.2,  we  currently  have  a  Java 
GUI  that  can  define  a  sparse  grid  with  a  given  total  physical  extent.  The  grid  for  the  screen  can  be 
iteratively  refined  and/or  extended  using  the  CPDF  technique.  Efforts  are  presently  underway  to 
recode  in  C++  and  integrate  the  technique  into  WaveTrain™  as  the  class  CPDFScreen  to  be  derived 


29 


from  the  base  class  wtPhaseScreen. 


(a)  (b) 

Figure  5.2:  The  Java-based  GUI  for  generating  phase  screens  based  on  the  CPDF  technique.  The 
course  resolution  screen  in  shown  in  (a)  and  the  fine  resolution  screen  in  shown  in  (b). 


5.2.2  Closed  Loop  Tracking  Performance 

A  crucial  aspect  of  the  numerical  simulation  of  optical  propagation  through  atmospheric  turbulence 
is  the  generation  of  random  phase  screens  with  the  correct  statistics.  The  most  popular  technique 
used  is  the  DFT  using  the  computationally  efficient  fast  Fourier  transform  (FFT)  algorithm.  It  is 
well  known  that  in  order  to  adequately  model  the  low  spatial  frequency  content  of  the  atmospheric 
disturbance,  one  must  either  generate  large  (much  larger  than  the  propagation  grid)  phase  screens  or 
“boost”  the  low  frequency  content  of  the  random  screens  (Zernike  boost  or  Sub-harmonics  are  two 
such  methods).  Often  times  an  argument  is  made  that  if  atmospheric  tracking  (tilt  compensation) 
is  implemented,  the  low  frequency  content  can  be  ignored  and  the  FFT  technique  is  adequate.  That 
argument  is  investigated  in  this  study.  Comparisons  of  relevant  statistics  are  shown  under  various 
tracking  conditions  (including  anisoplanatism)  using  random  phase  screens  known  to  be  be  lacking 
in  low  spatial  frequency  content  (FFT)  and  random  screens  with  the  proper  low  spatial  frequency 
content. 

Using  the  modified  FS  approach  described  in  Sec.  5.2.1,  we  desired  to  determine  the  point  at 
which  the  low  spatial  frequency  disturbances  become  important.  Under  the  frozen  flow  hypothesis, 
the  temporal  tilt  (Zernike  2  and  3)  PSD  can  be  expressed  in  terms  of  the  spatial  phase  PSD  and 
the  velocity  [1] 
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where  K  is  normalized  spatial  frequency,  and  /i  =  fD/2v,  is  normalized  temporal  frequency.  For 
a  given  temporal  frequency,  normalized  spatial  frequencies  below  2nfi  do  not  contribute  to  the 
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temporal  PSD.  When  using  phase  screens  with  a  minimum  spatial  frequency  content  of  K^niu  > 
2-7t/i,  normalized  temporal  frequencies  below  ^lrlin  are  under-represented.  An  example  temporal 
phase  spectrum  (from  a  wave-optics  simulation)  is  shown  in  Fig.  5.3  for  phase  screens  generated 
using  the  FFT  technique  and  the  modified  FS  approach.  As  can  be  seen  in  the  plot,  at  low  temporal 


Figure  5.3:  Example  temporal  spectra  for  two  different  wave-optics  simulations  each  using  a  different 
technique  for  generating  random  phase  disturbances. 


frequencies,  the  power  is  low  when  compared  to  the  theoretical  PSD. 

We  have  run  a  set  of  wave-optics  simulations  using  WaveTrain™  mex  systems  to  determine  the 
dependence  of  closed-loop  tracking  performance  on  the  spatial  frequency  content  of  the  phase  dis¬ 
turbances.  We  implemented  point  source  tracking  for  a  “tactical”  HEL  scenario  [see  Table  5.2].  The 


Table  5.2:  Simulation  Parameters  for  closed  loop  tracking  study 


Parameter 

Value 

Platform  Altitude 

3000  m 

Target  Altitude 

10  m 

Platform  Velocity 

100  m/s 

Target  Velocity 

10  m/s 

Aperture  Diameter 

0.5  m 

Wavelength 

1.064  (im 

Cl  Profile 

2xHV5/7 

Spherical  ro 

0.27  m 

Rytov 

0.11 

Isoplanatic  Angle 

1.5  /irads 

Greenwood  Frequency 

189  Hz 

Tyler  Frequency 

28  Hz 

Frame  Rate 

1000  Hz 

first  step  was  to  verify  the  error  rejection  function  of  our  tracker.  The  theoretical  error  rejection 
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function  for  a  simple  integrator  with  latency  is 


ERJ(/) 


1  + 


sin(27r/Af) 


-l 


(5.8) 


where  fsw  is  the  3dB  closed  loop  bandwidth,  and  At  is  the  latency, 
and  closed  loop  gain,  f3,  are  related  by 
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2n 


fsw 
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The  closed  loop  bandwidth 


(5.9) 


where  Fs  is  the  frame  rate.  We  set  up  a  simple  simulation  to  verify  this  form  of  the  error  rejection 
function.  The  WaveTrain™  system  and  the  results  of  the  error  rejection  test  are  shown  in  Fig.  5.4. 
As  can  be  seen  in  the  figure,  the  form  of  the  error  rejection  function  from  the  simulation  matches 


Disturbance  Frequency  (Hz)  Disturbance  Frequency  (Hz) 


(c) 


(d) 


Figure  5.4:  Verification  of  the  WaveTrain™  error  rejection  function,  (a)  The  WaveTrain™  system 
used  and  the  error  rejection  test  results  for  (a)  0  frames,  (b)  1  frame,  and  (c)  2  frames  of  latency. 
The  lines  are  the  theoretical  values  [see  (5.8)]  and  the  circles  are  the  results  from  the  WaveTrain™ 
simulation. 
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the  theoretical  form  from  (5.8). 

The  simulation  results  for  closed  loop  tracking  are  shown  in  Fig.  5.5  as  a  function  of  the  normal¬ 
ized  (by  Tyler  frequency)  closed  loop  bandwidth.  At  low  closed  loop  bandwidth  (less  than  the  Tyler 


Target  Board  -  0  Frames  Latency 


Normalized  BW 

Figure  5.5:  Target  board  residual  jitter  results  of  closed  loop  tracking  simulation  using  standard 
FFT  screens  (blue  circles)  and  modified  FS  screens  (green  circles)  as  compared  to  the  theoretical 
values  (black  line). 


frequency)  the  discrepancy  between  the  FFT  screens  and  the  FS  screens  is  obvious,  especially  for 
2  frames  of  latency.  However,  as  the  bandwidth  increases  to  greater  than  the  Tyler  frequency,  the 
resulting  residual  jitter  appears  to  be  independent  of  the  spatial  frequency  content  in  the  random 
screens.  In  order  to  get  more  than  90%  of  the  residual  jitter  predicted  by  the  FS  approach  using  a 
FFT  approach,  the  closed  loop  bandwidth  must  be  greater  than  the  Tyler  frequency. 


5.2.3  Orbital  Propagator 

A  recent  addition  to  the  Scaling  for  HEL  and  Relay  Engagements  (SHaRE)  Matlab®  toolbox  is  the 
capability  of  generating  engagement  scenarios  which  include  orbital  assets.  We  have  incorporated 
into  SHaRE  the  North  American  Aerospace  Defense  Command  (NORAD)  SGP4  (for  near-earth 
objects)  and  SDP4  (for  deep-space  objects)  algorithms  for  determining  satellite  location  and  velocity 
in  earth  orbit  using  current  NORAD  two-line  element  (TLE)  datum.  This  allows  us  to  investigate 
the  energy  characteristics  of  a  ground  based  laser  to  an  orbital  asset.  We  have  conducted  a  study 
using  such  a  simulation  to  demonstrate  the  use  of  the  orbital  propagator  under  a  diverse  set  of 
atmospheric  conditions.  The  ephemeris  data  for  selected  satellites  is  retrieved  from  freely  released 
TLE  satellite  catalogs  (www.space-track.org).  The  atmospheric  data  and  satellites  state  vectors 
are  then  feed  into  a  line-of-site  algorithm  for  analysis.  This  determines  the  access  times  and  look 
angles  between  ground  sites  and  assets.  Based  on  the  engagement  geometry,  SHaRE  was  used  to 
quantitatively  define  the  open  loop  and  closed  loop  characteristic  of  the  energy  at  target.  An  example 
of  computing  satellite  positions  and  visibility  from  a  ground  station  is  shown  in  Fig.  5.6. 
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Figure  5.6:  Satellite  ground  track  and  visibility  from  WPAFB  on  10  October  2006.  Solid  lines 
represent  times  at  which  the  satellite  would  have  been  visible  from  a  ground  station  at  WPAFB. 


5.2.4  Extreme  and  Percentile  Environmental  Reference  Tables  Atmo¬ 
sphere 


As  a  by-product  of  an  effort  to  integrate  the  SHaRE  Matlab®  toolbox  with  the  Air  Force  Institute 
of  Technology  (AFIT)  High  Energy  Laser  End-to-End  Simulation  (HELEEOS),  we  have  included 
AFIT’s  extensive  Extreme  and  Percentile  Environmental  Reference  Tables  (ExPERT)  based  atmo¬ 
spheric  parameter  data  into  SHaRE.  The  ExPERT  database  is  a  joint  effort  by  the  Air  Force 
Research  Laboratory’s  Air  Vehicles  and  Space  Vehicles  Directorates,  and  the  Air  Force  Combat 
Climatology  Center.  ExPERT  is  a  Microsoft  Access  database  of  pre-calculated  climatological  values 
for  various  regions-land,  ocean,  and  upper  air-as  well  as  408  sites  worldwide.  For  the  individual 
surface  land  sites,  ExPERT  allows  the  analyst  to  view  monthly  and  hourly  percentile  data,  duration 
data,  and  yearly  minimum  and  maximum  values  for  the  following  atmospheric  variables:  altimeter 
setting,  dew  point  temperature,  absolute  humidity,  relative  humidity,  specific  humidity,  tempera¬ 
ture,  wind  speed,  and  wind  speed  with  gusts.  Percentiles  for  diurnal  data  and  sky  cover  data  are 
displayed  as  well.  Also  available  are  the  percent  frequency  of  occurrence  for  several  “significant” 
weather  phenomena:  thunderstorms,  fog,  blowing  snow  or  sand,  freezing  rain,  hail,  snow,  and  rain. 
Notably,  ExPERT  also  enables  the  analyst  to  display  the  probabilities  of  when  a  particular  combi¬ 
nation  of  temperature  and  relative  humidity  will  occur  for  a  specific  land  site.  Using  ExPERT  and 
the  Global  Aerosol  Data  Set  (GADS),  along  with  the  High- Resolution  Transmission  Molecular  Ab¬ 
sorption  Database  (HITRAN),  AFIT  has  generated  a  look-up  table  (LUT)  for  molecular  and  aerosol 
absorption  and  scattering,  temperature,  and  pressure.  The  free  parameters  in  the  LUT  available 
for  use  in  ATMTools  are  altitude  (up  to  80  kft),  boundary  layer  altitude,  wavelength  (24  values), 
site  (408  worldwide),  season  (summer  or  winter),  time  of  day  (8  periods  and  a  daily  average),  and 
relative  humidity  (9  values).  Fig.  5.7  shows  the  available  locations  in  the  ExPERT  database. 
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Figure  5.7:  Available  locations  in  the  ExPERT  database.  The  sites  are  classified  into  5  categories, 
desert  (blues),  north  latitudes  (cyan),  polar  (red),  south  latitudes  (magenta),  and  tropical  (green). 


5.2.5  Implementation  of  Adaptive  Controller  in  WaveTrain™ 

The  adaptive  controller  developed  by  UCLA  was  implemented  in  Matlab®/ Simul  ink®  models  using 
mex  functions.  In  an  effort  to  incorporate  the  adaptive  controller  in  a  WaveTrain™  simulation,  MZA 
worked  with  UCLA  graduate  assistant  Yutai  Liu  to  first  call  the  Matlab®  functions  from  a  Wave- 
Train™  block  via  the  Matlab®  engine  (called  an  m-system)  and  then  later  to  have  the  WaveTrain™ 
block  directly  call  the  functions  in  the  C++  code  underlying  the  mex  functions.  The  implementation 
was  verified  by  comparing  the  WaveTrain™  and  Simulink®  results. 

5.2.6  Laboratory  and  Wave  Optics  Analysis  of  Aero- Optical  Effects 

Aero-optical  OPD  wind  tunnel  data  taken  at  the  University  of  Notre  Dame  (UND)  was  applied  to 
a  liquid  crystal  (LC)  spatial  light  modulator  (SLM)  at  the  AFIT  adaptive  optics  (AO)  laboratory 
at  Wright-Patterson  Air  Force  Base  (WPAFB).  Captain  Jason  Schmidt  (AFIT/ENG)  permitted 
MZA  to  use  this  facility  during  the  assembly  of  the  MZA  laboratory.  The  wind  tunnel  data  was  also 
incorporated  as  an  OPD  in  a  WaveTrain™  simulation.  Jitter,  Strehl,  and  wavefront  sensor  (WFS) 
slopes  were  analyzed.  Comparisons  showed  that  laboratory  and  WaveTrain™  results  agreed  well. 

Aero-optic  OPD  Data 

The  data  provided  by  the  UND  Aero-Optics  Group  included  two  runs  each  with  the  turret  at 
120°  and  130°  azimuth  angles.  Zernike  decomposition  of  the  OPD  data  on  a  frame-by-frame  basis 
was  completed  for  Zernike  modes  1  to  55  (using  Noll’s  numbering  convention  [2]).  The  PSD  was 
calculated  for  each  Zernike  mode.  Fractional  power  of  each  mode  (up  to  Zernike  10)  for  each  of  the 
four  runs  is  presented  in  Fig.  5.8.  Behavior  with  respect  to  disturbance  by  Zernike  mode  cannot  be 
generalized  other  than  noting  that  higher  order  Zernike  modes  have  more  fractional  power  in  the 
130°  data  sets  than  the  120°  data  sets.  PSDs  for  modes  with  total  power  greater  than  5%  of  the 
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total  disturbance  power  are  presented  in  Fig.  5.9.  For  both  the  120°  and  130°  data  the  disturbance 
power  is  dominant  in  frequencies  less  than  or  equal  to  100  Hz.  In  the  120°  data,  the  drop  in  power 
for  frequencies  greater  than  100  Hz  is  sharper  than  that  of  the  130°  data,  and  the  variation  in 
disturbance  frequency  modes  for  higher  frequencies  is  more  pronounced  in  the  130°  data. 


Figure  5.8:  Fraction  of  total  power  in  the  first  10  Zernike  modes  for  each  of  the  data  sets. 


WaveTrain™  Modeling 

A  WaveTrain™  model  for  simulating  the  effects  of  propagation  thorough  an  aero-optics  disturbance 
was  developed.  The  WaveTrain™  wave-optics  model  is  a  simple  model  which  included  the  OPD 
input  as  its  sole  disturbance.  Modeling  parameters  are  shown  in  Table  5.3.  Strehl  values  were 
calculated  using  a  Gaussian  quadrature  analysis  of  intensity  degradation.  Jitter  was  determined  as 
the  time-variance  of  the  centroid  positions. 


Table  5.3:  Modeling  parameters  for  WaveTrain™  aero-optics  model 


Parameter 

Value 

Telescope  Diameter 
Source  Wavelength 

7.55  mm 
632.8  nm 

Laboratory  Analysis 

For  the  experimental  portion,  we  used  AFIT’s  AO  laboratory.  The  setup  included  a  Helium:Neon 
(HeNe)  laser  at  632.8  nm,  a  512x512  SLM,  a  32  subaperture  Shack-Hartmann  WFS,  and  a  far-field 
scoring  camera.  The  illuminated  portion  of  the  SLM  and  WFS  measured  7.55  /rm  across  and  the 
camera  pixel  spacing  was  6.7  /rm.  The  OPDs  were  up-sampled  to  a  512x512  grid,  applied  to  the 
SLM,  and  corresponding  WFS  slopes  and  far-field  camera  images  were  recorded.  For  comparison,  a 
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Figure  5.9:  PSDs  for  Zernike  modes  which  represented  more  than  5%  of  total  disturbance  power  by 
data  file  (a)  120°,  Run  1,  (b)  120°,  Run  2,  (c)  130°,  Run  1  and  (d)  130°,  Run  2 


set  of  data  was  taken  with  a  zero  OPD  applied  to  the  SLM.  The  far-held  images  were  processed  by 
time-averaging  and  tilt-removed  time  averaging.  The  higher-order  Strehl  was  calculated  as  the  ratio 
of  the  peak  of  the  tilt-removed,  time-averaged  image  to  the  peak  value  of  the  comparison  zero  OPD 
set.  The  long-term  Strehl  values  were  calculated  as  the  ratio  of  the  peak  of  the  time-averaged  image 
to  the  peak  value  of  the  comparison  set.  Jitter  was  computed  by  taking  the  standard  deviation  of  the 
centroid  shift  on  a  frame-by-frame  basis.  The  laboratory  parameters  are  summarized  in  Table  5.4. 
The  calibration  run  with  zero  OPD  data  indicated  a  slight  change  in  SLM  behavior  over  time  of 


Table  5.4:  Laboratory  Parameters 


Parameter 

Value 

Telescope  Diameter 
Source  Wavelength 

0.28  m 

1  /im 

use.  The  centroid  drifts  over  a  sequence  of  450  frames.  This  indicates  an  effect  of  either  heating  or 
continuous  voltage  application.  This  centroid  drift  can  be  seen  in  Fig.  5.10. 
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Figure  5.10:  Indication  of  drift  for  SLM  use  over  time  for  zero  OPD  run 


Wavefront  Sensor  Comparison  with  Spatial  Light  Modulator  Input 

The  wavefront  was  reconstructed  from  the  WFS  slopes  and  compared  to  the  OPDs  placed  on  the 
SLM.  The  mean-squared  error  (MSE)  between  the  input  OPD  and  the  reconstructed  phase  was 
calculated  on  a  frame-by-frame  basis.  The  reconstructed  phase  from  the  WFS  agrees  very  well  with 
the  OPD  data  when  the  phase  variance  is  small  (120°,  Run  2).  The  agreement  is  less  accurate  for 
cases  with  greater  disturbances,  but  the  WFS  phase  variance  results  reflect  the  overall  trend  in  the 
OPD  variance  [see  Fig.  5.11].  Examples  of  the  best  match  (in  terms  of  MSE)  and  worst  match  phase 
reconstructions  and  time-average  far  field  irradiance  patterns  are  shown  in  Fig.  5.12  and  Fig.  5.13. 
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Figure  5.11:  Mean  squared  error  between  the  input  OPD  phase  and  the  reconstructed  phase  from 
WFS  slope  measurements  in  the  laboratory. 
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Figure  5.12:  Example  input  and  reconstructed  phase  maps  for  (a)  the  best  match  (frame  336  from 
120°,  Run2)  and  (b)  the  worst  match  (frame  23  from  130°,  Runl).  Even  in  the  worst  case,  some  of 
the  same  trends  can  be  observed  in  the  phase  reconstruction  as  in  the  input  OPD. 
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Figure  5.13:  Example  WaveTrain™  and  laboratory  measured  time-averaged  far  field  patterns  for  (a) 
120° ,  Run2  and  (b)  130°,  Runl. 


Results  Summary 

The  WaveTrain™  results  agree  very  well  with  the  laboratory  results  and  averaged  images  share 
distinctive  features  [see  Fig.  5.13].  Strehl  values  show  the  same  trends  and  very  similar  values.  The 
jitter  values  are  slightly  higher  in  the  laboratory  than  in  the  WaveTrain™  model.  This  difference  can 
be  partially  explained  by  the  SLM  drift  observed  with  the  zero  OPD  runs.  A  results  summary  is 
shown  in  Table  5.5. 
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Table  5.5:  Results  comparison  summary 


Data  Set 

120°,  Run  1 

120°,  Run  2 

130°,  Run  1 

130°,  Run  2 

WaveTrain™  Strehl 

0.21 

0.59 

0.31 

0.24 

Laboratory  L.T.  Strehl 

0.24 

0.60 

0.26 

0.19 

Laboratory  H.O.  Strehl 

0.25 

0.61 

0.25 

0.19 

WaveTrain™  X  Jitter  (A/D) 

0.0924 

0.0480 

0.0420 

0.0238 

WaveTrain™  Y  Jitter  (A/D) 

0.1512 

0.0258 

0.0219 

0.0177 

Laboratory  X  Jitter  (A/D) 

0.1272 

0.0904 

0.1777 

0.1312 

Laboratory  Y  Jitter  (A/D) 

0.1706 

0.1440 

0.1486 

0.2091 
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Chapter  5 

Active  Contours  for  Tracking 

Georgia  Tech 

5.1  Introduction 

An  active  contour  approach  driven  by  ideas  from  knowledge-based  segmentation  was  developed  for  missile 
tracking  through  deep  turbulence.  Prior  work  on  laser-guided  target  tracking  for  missiles  has  focused  on 
understanding  the  source  of  jitter  and  reducing  jitter  through  controller  design  for  the  optical  system. 
Although  effective  plant  and  controller  design  lead  to  performance  improvements,  turbulence  effects  are  a 
consistent  source  of  noise  arising  from  the  image  processing  not  modelled  within  the  physics  underlying  tilt 
stabilization,  which  is  what  an  effective  tilt  controller  design  typically  corrects  for. 

Our  research  has  shown  how  geometric  active  contours  driven  by  Bayesian  statistical  techniques  may  be 
used  in  order  to  dynamically  track  missiles  through  deep  turbulence. 

We  can  summarize  our  contributions  as  follows: 

•  All  of  the  smoothing  we  do  for  the  posteriors  is  based  on  active  contour  methods  (anisotropic  image 
enhancement).  Thus  each  isophote  (equal  probability  contour)  is  moving  according  to  a  standard  curve 
evolution  flow. 

•  Active  contours  without  edges  (region-based  active  contours)  appear  in  the  estimation  of  the  gain 
factor  for  separating  the  foreground  from  the  background  classes. 

•  Active  contours  are  used  for  tip  tracking. 

We  now  briefly  summarize  some  key  aspects  of  our  work  in  the  next  sections. 


5.2  Bayesian  Tracking 

The  active  contours  are  driven  by  a  Bayesian  estimation  technique  as  follows:  Suppose  that  each  thresh¬ 
old  graylevel  belongs  to  a  particular  class,  the  image  to  segment  consists  of  M  classes  of  graylevels  C  = 
{1 ,M},  and  that  each  class  c  £  C  is  associated  with  a  mean  pixel  intensity  \xc  and  standard  deviation 
trc.  The  goal  of  the  segmentation  process  is  to  map  each  pixel  intensity  into  a  class.  A  pixel  intensity 
measurement  17  is  a  random  variable  with  known  distribution,  independent  of  the  other  pixels  in  the  image 
domain.  Given  prior  knowledge  of  the  image  statistics,  the  probability  of  a  pixel  intensity  being  assigned  to 
a  particular  class  is  given  by  Bayes’  Rule, 


Pr(ci  =  c\vi  =  v ) 


Pr(uj  =  v\cj  =  c)  Pr(cj  =  c) 

Z)7  Pr(ty  =  v\ci  =  7)  Pr(ci  =  7)  ’ 


(5.1) 
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where  the  probability  Pr(i>i  =  i>|cj  =  c)  is  assumed  to  have  a  Gaussian  distribution1, 

Prfa  =  v\a  =  c)  =  Jl  exp{— — ^-},  (5.2) 

V27T  <rc  2  crj 

with  (//c,  crc)  the  mean  and  standard  deviation  corresponding  to  the  class  c  £  C.  The  posterior  probabilities 
from  equation  (4.1)  are  calculated  for  each  class,  then  smoothed  and  normalized  to  obtain  Pr.  The  pixel  i 
is  assigned  to  the  class  with  the  maximal  smoothed  posterior  probability, 

c*  =  argmaxPr(cj  =  c\vi  =  v).  (5-3) 

cgC 

Computing  the  posterior  probabilities  begins  with  the  assumption  of  known  image  statistics,  i.e. ,  the 
prior  probabilities  Pr(c.j  =  c)  are  needed,  as  are  the  mean  and  standard  deviations  corresponding  to  the 
classes  C.  Initially  the  probabilities  Pr(cj  =  c)  are  assumed  to  be  homogeneous.  After  the  first  image,  the 
priors  are  set  to  be  the  posteriors  from  the  previous  time  step.  The  means  /. ic  and  standard  deviations  ac  for 
the  classes  c  £  C  are  updated  from  their  previous  state  according  to  the  computed  statistics  of  the  current 
segmentation. 

Implementing  Bayes’  Rule  for  a  sequence  of  images  without  smoothing  of  the  posteriors  will  tend  to  result 
in  convergence  to  a  fixed  steady  state  segmentation,  potentially  incorrect.  Such  a  convergence  is  undesirable 
when  dealing  with  moving  imagery,  as  is  the  case  discussed  here.  To  prevent  the  segmentation  process  from 
converging,  smoothing  of  the  posteriors  is  performed.  The  smoothing  process,  performed  prior  to  the  MAP 
classification,  has  the  additional  benefit  of  removing  noise. 


5.3  Background  on  Curve  Evolution 

One  can  think  of  an  image  as  a  map  /  :  2)  — »  £,  i.e.,  to  any  point  x  in  the  domain  2),  /  associates  a 
“color”  I(x)  in  a  color  space  £.  For  ease  of  presentation  we  will  mainly  restrict  ourselves  to  the  case  of  a 
two-dimensional  gray  scale  image  which  we  can  think  of  as  a  function  from  a  domain  2)  =  [0,  1]  x  [0, 1]  C  R2 
to  the  unit  interval  £  =  [0, 1]. 

The  algorithms  all  involve  solving  the  initial  value  problem  for  some  PDE  for  a  given  amount  of  time. 
The  solution  to  this  PDE  can  be  either  the  image  itself  at  different  stages  of  modification,  or  some  other 
object  (such  as  a  closed  curve  delineating  object  boundaries)  whose  evolution  is  driven  by  the  image. 

For  example,  introducing  an  artificial  time  t,  the  image  can  be  deformed  according  to 

%  =  ™  <5-4> 

where  I(x,t)  :  2)  x  [0,  T)  — >  £  is  the  evolving  image,  T  is  an  operator  which  characterizes  the  given  algorithm, 
and  the  initial  condition  is  the  input  image  Iq.  The  processed  image  is  the  solution  I(x,t)  of  the  differential 
equation  at  time  t.  The  operator  T  usually  is  a  differential  operator,  although  its  dependence  on  I  may  also 
be  nonlocal. 

Similarly,  one  can  evolve  a  closed  curve  r  C  2)  representing  the  boundaries  of  some  planar  shape  (r  need 
not  be  connected  and  could  have  several  components).  In  this  case,  the  operator  T  specifies  the  normal 
velocity  of  the  curve  that  it  deforms.  In  many  cases  this  normal  velocity  is  a  function  of  the  curvature  k  of 
r,  and  of  the  image  /  evaluated  on  T.  A  flow  of  the  form 

^=^(7,«)N  (5.5) 

is  obtained,  where  N  is  the  unit  normal  to  the  curve  T. 

1  Alternative  distributions  better  suited  to  the  image  statistics  may  be  picked  for  the  classes.  A  more  representative  distri¬ 
bution  improves  the  segmentation  process. 
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Very  often,  the  deformation  is  obtained  as  the  steepest  descent  for  some  energy  functional.  For  example, 
the  energy 

m  =IJ  II  V/f  dxdy  (5.6) 

and  its  associated  steepest  descent,  the  heat  equation, 


dl 

dt 


=  A I 


(5.7) 


correspond  to  the  classical  Gaussian  smoothing. 

The  use  of  PDEs  allows  for  the  modelling  of  the  crucial  but  poorly  understood  interactions  between 
top-down  and  bottom-up  vision.  In  a  variational  framework,  for  example,  an  energy  £  is  defined  globally 
while  the  corresponding  operator  T  will  influence  the  image  locally.  Algorithms  defined  in  terms  of  PDEs 
treat  images  as  continuous  rather  than  discrete  objects.  This  simplifies  the  formalism,  which  becomes  grid 
independent.  On  the  other  hand  models  based  on  nonlinear  PDEs  may  be  much  harder  to  analyze  and 
implement  rigorously. 


5.4  Geometric  Active  Contours 

For  geometric  active  contours,  one  deforms  the  active  contour  T  by  a  velocity  which  is  essentially  defined 
by  a  curvature  term,  and  a  constant  inflationary  term  weighted  by  a  stopping  function  W.  By  formulating 
everything  in  terms  of  quantities  which  are  invariant  under  reparametrization  (such  as  the  curvature  and 
normal  velocity  of  T)  one  obtains  an  algorithm  which  does  not  depend  on  the  parametrization  of  the  contour. 
In  particular,  it  can  be  implemented  using  level  sets. 

More  specifically,  a  simple  geometric  model  is  given  by 

V  =  W(x)(k  +  c),  (5.8) 

where  both  the  velocity  V  and  the  curvature  k  are  measured  using  the  inward  normal  N  for  T.  Here,  as 
previously,  W  is  small  at  edges  and  large  everywhere  else,  and  c  is  a  constant,  called  the  inflationary  param¬ 
eter.  When  c  is  positive,  it  helps  push  the  contour  through  concavities,  and  can  speed  up  the  segmentation 
process.  When  it  is  negative,  it  allows  expanding  “bubbles,”  i.e.,  contours  which  expand  rather  than  contract 
to  the  desired  boundaries.  We  should  note  that  there  is  no  canonical  choice  for  the  constant  c,  which  has  to 
be  determined  experimentally. 

In  practice  T  is  deformed  using  the  Osher-Sethian  level  set  method.  One  represents  the  curve  Tt  as  the 
zero  level  set  of  a  function  <f>  :  £>  x  M+  — >  R, 


r‘  =  {xG  D  :  $(x,f)  =  0}.  (5.9) 

For  a  given  normal  velocity  field,  the  defining  function  <I>  is  then  the  solution  of  the  Hamilton- Jacobi 
equation 

9$  M  .. 

— +  V||V0>||=0 

which  can  be  analyzed  using  viscosity  theory. 

Geometric  active  contours  have  the  advantage  that  they  allow  for  topological  changes  (splitting  and 
merging)  of  the  active  contour  T.  The  main  problem  with  this  model  is  that  the  desired  edges  are  not 
steady  states  for  the  flow  4.8.  The  effect  of  the  factor  W(x)  is  merely  to  slow  the  evolution  of  Tt  down  as  it 
approaches  an  edge,  but  it  is  not  the  case  that  the  Tt  will  eventually  converge  to  anything  like  the  sought-for 
edge  as  t  — ■>  oo.  Some  kind  of  artificial  intervention  is  required  to  stop  the  evolution  when  rt  is  close  to  an 
edge. 
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5.4.1  Conformal  (Geodesic)  Active  Contours 

For  conformal  active  contours,  one  defines  a  Riemannian  metric  gw  on  £)  from  a  given  image  I  :  3  — >  K,  by 
conformally  changing  the  standard  Euclidean  metric  to, 

gw  =  W{x)2\\dx\\2.  (5.10) 

The  length  of  a  curve  in  this  metric  is 

CW{T)  =  J  W(r(s))  ds.  (5.11) 


Curves  which  minimize  this  length  will  prefer  to  be  in  regions  where  W  is  small,  which  is  exactly  where  one 
would  expect  to  find  the  edges.  So,  to  find  edges,  one  should  minimize  the  IT-weighted  length  of  a  closed 
curve  r,  rather  than  some  “energy”  of  T  (which  depends  on  a  parametrization  of  the  curve). 

To  minimize  £lv(T),  one  computes  a  gradient  flow  in  the  L2  sense.  Since  the  first  variation  of  this  length 
functional  is  given  by 


d  CW{Y) 
At 


N  •  VIE}  ds, 


where  V  is  the  normal  velocity  measured  in  the  Euclidean  metric,  and  N  is  the  Euclidean  unit  normal,  the 
corresponding  L 2  gradient  flow  is 

Uconf  =  Wk  —  N  •  VW.  (5.12) 


Contemplation  of  the  conformal  active  contours  leads  to  another  interpretation  of  the  concept  “edge.” 
Using  the  landscape  metaphor  one  can  describe  the  graph  of  IT  as  a  plateau  (where  | 
nablal |  is  small)  in  which  a  canyon  has  been  carved  (where  |V/|  is  large).  The  edge  is  to  be  found  at  the 
bottom  of  the  canyon.  Now  if  IT  is  a  Morse  function,  then  one  expects  the  “bottom  of  the  canyon”  to  consist 
of  local  minima  of  W  alternated  by  saddle  points.  The  saddle  points  are  connected  to  the  minima  by  their 
unstable  manifolds  for  the  gradient  flow  of  W  (the  ODE  af  =  —  X/W(x).)  Together  these  unstable  manifolds 
form  one  or  more  closed  curves  which  one  may  regard  as  the  edges  which  are  to  be  found. 


5.5  Region-Based  Active  Contours 

We  can  easily  include  statistic  globally  based  terms  in  our  framework.  This  is  derived  from  the  following 
type  problem: 

minimize  /  fidx  +  /  f2da:  (5.13) 

./interior  C  ./exterior  C 

for  globally  defined  functions  fi  and  /2.  For  example,  to  separate  means,  we  define  u  =  Su/Au,  w  =  Sw/Aw, 
and 


Su  = 

[  I  cla: 

Au  — 

f  Ax 

Jr " 

Jr “ 

sw  = 

[  I  Ax 

Aw= 

/  dx, 

Jr ™ 

lRw 

and  Ru  and  Rw  denote  the  domains  inside  and  outside  the  curve  respectively.  Computing  the  first  variation 
results  in  the  following  gradient  flow  for  separating  the  means  of  regions 

rt  =  {u-w)  (I—r^  +  ~T—)  N-  (5-14) 
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